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INTRODUCTION 


A  substantial  need  has  been  demonstrated  within  the  Armament  Research, 
Development  and  Engineering  Center  (ARDEC)  for  a  convenient-to-use  and 
computationally  efficient  flight  simulation  for  smart  munitions.  To  meet  this 
need,  ARDEC  is  formulating  methods  and  writing  modular  software  for 
rapidly  developing  autopilot  simulations  for  guided  projectiles  and  other 
munitions. 

Time  constants  associated  with  autopilot  components  are  often  small 
compared  with  their  driving  terms.  The  integration  time  step  is  consequently 
driven  to  very  small  values  to  achieve  stable  numerical  integration,  which 
results  in  increased  computer  run  time.  An  innovative  technique  was 
developed  in  which  exact  analytic  solutions  to  differential  equations  and 
transfer  functions  are  applied  in  a  piecewise  manner  within  a  larger  but  lower 
frequency  problem  that  is  solved  numerically.1 

Closed  form  analytical  solutions  were  previously  obtained  for  the 
following:  the  first-order  lag,  the  first-order  lag  with  differentiation,  the 
first-order  lead/lag,  the  first-order  lag  with  integrator,  the  second-order 
lag/oscillator,1  a  two-axis  gimbaled  gyro,2  and  an  impulse  thruster.3  The 
typical  autopilot  can  be  described  as  a  concatenation  of  these  functions  and 
rate  sensors,  switches,  limiters,  and  dead  zones.  When  small  integration  time 
steps  are  used,  as  in  the  numerical  integration  of  transfer  functions,  it  is 
straightforward  to  concatenate  solutions  in  serial  order  since  the  sampling  rate 
of  the  signal  is  high  because  of  the  small  integration  time  step.  However,  if 
analytical  solutions  to  transfer  functions  are  to  be  used  to  make  possible  the 
use  of  relatively  large  time  steps,  the  sampling  rate  may  no  longer  be 
adequate  to  sample  the  frequency  content  of  the  signal  as  it  moves  from  one 
transfer  function  to  another. 

This  report  discusses  the  replacement  of  serial  configurations  of  transfer 
functions  by  equivalent  parallel  configurations.  This  approach  avoids 
difficulties  arising  from  propagation  of  the  signal  through  a  sequential 
network  of  widely  varying  natural  frequencies  when  using  a  relatively  large 
piecewise  integration  time  step,  and  produces  a  decomposition  that  always 
leads  to  terms  that  can  be  readily  integrated  analytically. 


1  Eugene  M.  Friedman  and  Michael  J.  Amoruso,  "An  Analytical  Modularized  Treatment  of  Autopilots  for  Guided 
Projectile  Simulations,"  US  Army  Armament  Research  and  Development  Center,  Technical  Report  ARLCD-TR- 
85025,  Picatinny,  NJ,  August  1985. 

2 

M.  J.  Amoruso,  "An  Analytical  Technique  for  Modeling  Gyroscopes  in  Guided  Projectile  Simulations,"  US  Army 
Armament  Research  and  Development  Center,  Technical  Report  ARAED-TR -86010,  Picatinny,  NJ,  June  1986. 

3 

Eugene  M.  Friedman,  Michael  J.  Amoruso,  and  Romel  Campbell,  "An  Analytical  Model  of  an  Impulsive 
Thruster,"  Technical  Report  ARAED-TR-86034,  Picatinny  Arsenal,  NJ,  November  1986. 
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DISCUSSION 


When  modeling  guided  projectiles  in  6  DOF  (six  degree  of  freedom) 
simulations,  differential  equations  are  obtained  that  describe  the  various 
component  subsystems.  These  are  then  typically  integrated  numerically 
within  the  framework  of  the  6  DOF  simulation.  The  largest  allowable  time 
step  to  perform  the  integration  is  bounded  by  two  constraints.  The  driving 
term  or  input  must  not  vary  appreciably  during  a  time  step  and  the  time  step 
must  be  sufficiently  small  to  insure  a  stable  integration. 

Since  the  driving  term  rates  are  commensurate  with  the  airframe  motion 
rates,  inherently  slower  processes  than  those  associated  with  the  autopilot, 
stable  integration  is  a  lower  bound  to  integration  step  size  than  the 
requirement  for  driving  terms  that  remain  essentially  constant  during  the 
integration  time  step.  By  using  analytic  closed-form  solutions  for  the 
differential  equations  for  the  autopilot,  the  second  constraint  appears  to  be 
eliminated.  An  innovative  technique  was  developed1,4,5  by  which  exact 
analytical  solutions  to  the  required  transfer  functions  are  applied  in  a 
piecewise  manner  within  a  larger  but  lower  frequency  problem  which  must  be 
solved  numerically.  The  use  of  these  piecewise  analytical  solutions  to  the 
transfer  functions  guarantees  valid  integration  of  the  autopilot  transfer 
functions  regardless  of  the  integration  time  step  (app  A). 

The  overall  system  is  usually  analyzed  into  simpler  terms  fn  sequential 
order  that  are  separately  solved  by  numerical  integration  techniques, 
assuming  that  the  input  or  driving  term  is  essentially  constant  during  the 
integration  time  step.  These  factors  are  concatenated  with  the  output  of  one 
factor  or  block  becoming  the  input  to  the  next  block.  Since  the  integration 
time  step  is  generally  quite  small  to  achieve  stable  numerical  integration, 
negligible  errors  are  introduced  as  the  signal  propagates  through  the  usually 
modest  number  of  transfer  function  blocks  down  to  the  output. 

This  approach  was  modified  by  introducing  analytical  closed-form 
solutions  for  the  transfer  function  factors  that  were  previously  treated 
numerically.  This  approach  permitted  large  time  steps  across  the  individual 
transfer  functions  and  yet  guaranteed  stable  intergration  contingent  only  upon 
the  input  remaining  essentially  constant  during  the  integration  time  step. 


4  M.  J.  Amoruso,  T.  F.  DeYoung,  D.  F.  Ladd,  and  R.  D.  Schulz,  "A  Comprehensive  Digital  Flight  Simulation  of 
the  Cannon-Launched  Guided  Projectile,”  R-TR-77-007,  U  S.  Army  Armament  Command,  Rock  Island,  January 
1977. 

5  M.  J.  Amoruso  and  D.  D.  Ladd,  "TELUM,  A  Comprehensive  Digital  Flight  Simulation  of  the  COPPERHEAD 
Projectile,”  ARLCD-SP-82003,  Picatinny  Arsenal,  NJ,  June  1982. 
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However,  if  several  factors  are  concatenated  to  represent  a  more  complex 
transfer  function,  the  final  output  can  be  found  to  depend  on  the  order  of  the 
transfer  function  factors.  This  difficulty  arises  because,  although  the  input  of 
a  block  might  be  essentially  constant  during  an  integration  time  step,  the 
output  might  not  be  if  the  frequency  response  of  the  block  is  relatively  high. 
This  output  becomes  the  input  to  the  next  transfer  function  block.  The 
requirement  that  the  input  to  this  next  block  be  essentially  constant  can  break 
down  unless  the  sampling  rate  is  high.  This  requires  a  smaller  integration 
time  step  and  longer  computer  execution  time.  This  dilemma  does  not  arise 
in  the  former  numerical  integration  approach  because  the  very  fine  time  step 
required  avoided  difficulties  associated  with  incompatibilities  of  bandwidth 
and  frequency  content  as  the  signal  propagated  from  block  to  block. 

The  approach  finally  used  was  the  conversion  of  a  complex  transfer 
function  to  a  parallel  representation  instead  of  a  serial  representation.  This 
technique  not  only  saves  computer  run  time  but  can  make  the  task  of 
programming  simpler  and  nearly  automatic. 

The  obvious  advantage  to  an  equivalent  parallel  representation  is  that  each 
block  receives  the  same  input  at  the  same  time  and  each  block  produces  its 
output  at  the  end  of  the  same  single  time  step.  These  outputs  do  not  become 
inputs  to  other  autopilot  transfer  function  blocks,  but  are  instead  summed  to 
produce  the  overall  output  for  the  overall  transfer  function.  Generally,  this 
technique  not  only  avoids  simulation  errors  arising  from  time  step  size 
incompatibilities  with  bandwidth,  internal  lag,  and  frequency  content  of  the 
signal  propagating  through  the  sequence  of  transfer  function  blocks  but  also 
has  additional  benefits.  Noise  should  be  treated  in  a  manner  that  is  consistent 
with  the  analytical  treatment  of  the  signals.  The  treatment  of  noise  can  be 
considerably  simplified  by  using  a  parallel  representation  of  the  transfer 
function,  for  reasons  analogous  to  those  adduced  above.  In  addition  to 
producing  an  algorithm  that  is  easy  to  apply  and  considerably  faster  than 
previous  numerical  integration  approaches,  parallel  decomposition  generally 
leads  to  a  combination  of  elementary  expressions,  whose  Laplace  inverse  and 
analytical  integration  arc  well  known. 

A  specific  example  of  this  approach  is  given  below  to  clarify  this 
approach.  A  general  treatment  is  developed  in  the  next  section  and  computer 
software  to  implement  this  technique  is  to  be  found  in  the  appendixes. 

This  approach  consists  of  the  following  steps: 

(1)  Writing  down  the  Laplace  operator  expression  including 
non-zero  initial  conditions  from  the  differential  equation 
description  or  from  the  block  diagram  (which  usually  will  not 
show  the  initial  conditions) 

(2)  Factoring 

(3)  Making  a  forma!  partial  fraction  expansion 
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(4)  Finding  the  expansion  coefficients 

(5)  Writing  down  the  expanded  Laplace  transform. 


The  latter  can  then  be  inverted  from  standard  tables  of  inverse  Laplace 
transforms  to  obtain  the  analytic  solutions  in  parallel  decomposition  or 
calculated  using  the  residue  theorem  of  complex  variables. 

It  is  worthwhile  to  emphasize  the  first  of  the  steps  enumerated  above. 
Autopilots,  seekers,  control  actuators,  and  other  components  of  guided 
projectiles  are  conventionally  described  in  block  diagrams  in  terms  of  the 
Laplace  operator  s.  Typically,  these  block  diagrams  represent  the  underlying 
differential  equations  only  if  the  initial  conditions  vanish.  This  is  a 
convenient  shorthand  notation  but  can  also  be  a^ource  of  confusion.  Recall 

that  for  initial  conditions  v  =  v(f  =  0),  v  =  (t  =  0),  etc.,  the  Laplace 

O'  '  O 

transform  of  the  derivative  of  a  time  domain  function  is  given  by  6,7 


L[y(t )  J  =  s L(y(f )]  -  y0 

(1) 

L[y  (t )]  =  s2 L[y(t )]  -  sy0->:0 

(2) 

or,  in  general,  where  y ^  ^  represents  the  n  derivative  of  y(t)  with  respect 
to  t, 


(«) 


n-l 


L Lv  ;]  =  snL[y  (f )]  -  sn' 


n-2  • 

>VS  -vo 


••  -  y, 


(n  -  1) 


(3) 


This  will  be  illustrated  in  the  example  given  below.  There  are  two  methods 
for  implementing  this  procedure.  The  first  involves  the  taking  of  limits  and 
derivatives  and  requires  the  factoring  of  the  transfer  functions  into  first-order 
systems.  The  second  requires  the  solution  of  sets  of  simultaneous  equations 
and  does  not  involve  limits  or  derivatives.  Factoring  the  transfer  function 
into  first-order  terms  is  optional  in  the  second  method.  These  two 
approaches  are  illustrated  below  by  a  concrete  example. 


^  Erwin  Kreysig,  Advanced  Engineering  Mathematics  ,  Fourth  Edition,  John  Wiley  and  Sons,  New  York,  NY, 
1979. 

^  Equation  (1)  can  be  derived  as  follows.  Let  £  —  y  ~  y  .  Thus  £(  t  —  0)  =  =  0  and  £  =  y.  Then 

L(y)  =  L(£)  =  ,VL(0  —  SL(y)  —  Sh(y)  =  ,VL(y)  —  y  .  Recall  L(A)  =  A/s  if  A  is  a  constant. 
A  similar  derivation  holds  for  higher  orders. 


First  Method:  Limit  Approach 

Step  1  -  Obtain  Laplace  Transform:  Consider  an  autopilot  component 
represented  by  the  following  block  diagram  in  figure  1.  The  differential 
equation  corresponding  to  this  block  diagram  is 


y(f)  +  y(t)  —  2  y(r)  =  Driving  term 
=  KU-tQ)  +  KQ  =  Kt  +  L 


(4) 


-> 


1 


s  +  s  —  2 


-> 


Figure  1.  Kxample  of  block  diagram 


The  driving  term  will  be  taken  to  be  K  (t  - 1 +  K Q  =  Kt  +  L  in  this  example, 
wheje  K  and  L  are  constants  and  t  is  time.  This  has  the  Laplace  transform 
K Is  +  Us.  Let  L{ y( t) }  =  Y(s).  The  Laplace  transform  of  the  differential 
equation  is  not 

2  K  L 

j  F(s)  +  5K(.r)  -  2K(.v )  =  —  +  —  (5) 


but  is  rather 

-  .vy0  -  y0  +  v Y ( ,v )  -  y()  -  2 Y(s) 


or 

2  A'  +  sL 

i  s  +  s  -  2  )Y(s)  -  y  (  1  +.v  )  -  y  =  -  (7) 

2 

S 


K  L 

=  —  +  —  (6) 


where  the  additional  ^erms  arc  due  to  the  initial  conditions  defined  by 


>0=  Y  i  =  and  >'0  = 


.  This  may  be  written 


dt  !'=  L 
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Y(s)  = 


(8) 


K  +  sL  +  .V  (y0  +  y())  +  .?  y0 


2  2 

.9  (.9  +.v  —  2) 


Step  2  -  Factor:  This  may  be  factored 


K  +  sL  +  .9  (y0+v0)  +  .9  y0 


Y  (  .9  ) 


.9  (a+2)(.9—  1) 


(9) 


Step  3  -  Expansion:  This  may  be  formally  expanded  as  follows: 


ABC  D 

Y  ( .9  )  =  +  “ '  +  + 

-  S  .9  +  2  .9—1 


(10) 


Step  4  -  Find  Expansion  Coefficients:  The  partial  fraction  expansion 
coefficients  A  through  D  are  now  evaluated. 

hind  A: 


lim  .9  Y ( .9 )  =  lim 

.(  -0  ,!--0 


K  +  sL  +  .v2(y0  +  y0)  +  .9  3y  0 


(  .9  +  2  )  (  .9  ~  1  ) 


[  B  C 
A  +  lim  j  + 

s  •()  I  .9 


D 


+ 


.9  +  2  .9 


K 

2 


Find  B : 


£_  \  2  t 

lim  [  9  Y  is )  J  =  lim 
i  ■(>  d.s  5  .0  "  (Is 


K  +  sL  +  .92(y0  +  y0)  +  ,93y( 


(  .9  +  2  )  (  .9  -  1  ) 


K 

2 


(11) 


(12) 
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This  transition  from  parallel  decomposition  or  expansion  is  shown  in  figure  2 
for  this  simple  case.  It  is  a  simple  matter  to  invert  this  expression  into  the 
time  domain.  Tables  of  inverse  Laplace  transforms  can  be  found  in  standard 
math  tables.8  Here,  the  inversion  of  (15)  is  given  by 


Kt 

(K  ~r2L) 

4y(,-4y0—  K +2L  ' 

-  2t 

K  +  L  +2y0+y0 

2 

4 

12 

e  + 

3 

(16) 


Q 

Charles  D.  Hodgman,  ed.,  C.R.C.  Standard  Mathematical  Tables,  Chemical  Rubber  Publishing  Co.,  Cleveland, 
OH,  1952.  . 
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(a)  Step  1 


> 


K  +  sL  +  s2(y  +>-  )  +  s^y 


2  2 
X  (  X 


4-  x 


2  ) 


-> 


(b)  Step  2 


-> 


K  +  sL  +  s2(yQ  +  yQ )  +  .v3>0 


■v  (  *  +  2  )  (  .v  -  1  ) 


-> 


(c)  Steps  3-4 


4>-p—  *y0-K  +  2L 
12(.9  +  2) 


K  +  L  +  2y0+y()  I 
3(-f  ~  1)  " 

O - > 

i  =~k  i 


~(K  +  2L) 
4  s 


Figure  2.  Example  of  transition  to  parallel  representation 
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Second  Method:  Algebraic  Approach 

The  first  step  in  the  algebraic  approach  is  always  the  same  as  above.  The 
second  and  third  steps  may  be  done  as  before,  i.e.,  factoring  the  denominator 
of  the  transfer  function  into  monomials  and  expanding  in  terms  of  these 
monomials.  Later  higher  order  terms  will  be  retained. 

Step  4  -  Find  Expansion  Coefficients:  Write  the  right  side  of  (10)  with  a 
least  common  denominator  and  equate  to  (9). 

K  +  sL  +  .v2(y0  +  y0)  +  /y0 

Y(s)  = - 

2  2 

.9  (s  +  .9  —  2) 

ABC  D 

=  —  +  —  +  -  + - 

v2  .9  (.9  +  2)  .9-1 

A  (.9  +.9  —  2)  +  Bs(s  +.9  —  2)  +  C S  2( .9  1 )  +  D.92(.9  +  2) 

=  (17) 

2  2 

.9  (.9  +  .9  -  2) 

The  denominators  are  equal.  Therefore  the  numerators  are  also  equal. 
Collecting  powers  of  s  yields 

K  ~sL-s'(y0  +  yl))+siy0  =  s~'[B  +  C  +  D]  + s"[A  +  B  -  C  +2D]  + s[A -2B]-2A  (18) 


Since  the  powers  of  s  are  linearly  independent,  this  is  equivalent  to  the 
following  set  of  equations  which  must  be  solved  simultaneously. 


>'o 

= 

B 

+  C 

+ 

D 

(19) 

>'o 

+ 

=  A 

+ 

B  -  C 

+  2D 

(20) 

L 

A 

-  2  B 

(21) 

K 

- 

2  A 

(22) 

solution 

of 

th 

ese 

equations 

for  the  coefficients 

A  through  D  is 

elementary.  The  values  of  the  coefficients  are  identical  to  the  previous  result, 
equatu  ns  (11)  through  ( 14) . 

Alternatively,  the  investigator  may  wish  to  retain  some  higher  order  terms 
rather  than  reduce  all  the  denominators  to  monomials,  perhaps  to  retain  a 
physical  interpretation  of  terms.  Note  that  step  1  in  figure  2  contains  a 
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second  order  term,  viz.,  (5  +  5  —  2).  It  is  possible  to  make  a  partial  fraction 
expansion  without  going  to  a  decomposition  that  only  includes  first  order 
terms.  This  can  be  illustrated  by  the  following  example.  Replace  (10)  by  a 
more  general  expression,  viz. 


V(s) 

Y(s)  =  - 


G(s) 

1  C  n  P  +  sQ 

<  i  j 

=  2  +  2 

D.  E . 

i  =  l  «  ;  =  /-!  ] 


(23) 


where  the  D.  are  first  ord^r  terms  of  the  form  s  —  r  and  the  E.  are  second 
order  terms  of  the  form  (.r  +  015  + (3).  This  may  be  expanded 

V(s) 

Y (s )  =  -  (24) 

G(s) 


C  D  ..D  E  ..£  ~C  D  D  ..£  -..IP  -  sQ  1 D  ..D  E  .  .E  +  ..[£  ~  s0  ]D  D  E  .  E 

1  2  /1*1  «  2  13  n  1  /  ♦  1  -  l1  l  I  t  -  2  n  *’L  «  J  1  /  /  *  1  i,  -  1 


The  denominators  are  equal.  Therefore  this  expression  can  be  true  if  and 
only  if  the  numerators  are  equal.  For  the  example  given  in  this  report 

C1  C2  P3+^3 

)'(.!)  =  -  +  -  +  -  (25) 

2  2 
.v  s  +s-2 

Cj(s  +5-2)  +  C25(.v2+5-2)  +  P3s  +  s3Q3 
2  2 

S  (5  +  .V  —  2) 

equating  the  numerator  of  (25)  to  the  numerator  of  (9)  and  collecting  similar 
powers  of  s  gives 


K  +  sL  +  s  (y0+y0)  +  s  yQ  =  s  (Q3+C:)  +  s  (Cj+Cj+Pj)  +  s(Cl~2C2)  ~  2.C {  (26) 


Equating  similar  powers  of  s  gives  the  following  simultaneous  equations: 


>>o=  23+  C2 


(27) 
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(28) 


-V0  +  ^0  =  C1  +  C2  +  P3 

L  =  Cl  ~  2C2 

K  =  -2 C1 

When  these  are  solved 

K 

2 

(K  +  2L) 

C2=  ~  - 

4 

(3K  + 2L) 

p3  =  +  >'o  +  : 

4 

(4yo  +  K  +  2L) 

Q,= 

4 


(29) 

(30) 

(31) 

(32) 

(33) 

(34) 


Thus 

K  (K  +  2L)  3/C  +  2L +  4(y0  +  y0)  (4^0  +  K  +  2L  )s 

Y(s )  =  -  -  -  -  +  -  +  -  (35) 

2  j2  4(^2  +  5  — 2)  4(.f2+j  —  2) 

When  Laplace  inverted,  (35)  becomes 

Kt  (K+2L)  3K  +  2L  +  4(y0+);o)  ( 4y0+K+2L ) 

y(l)  -  —  —  +  [e  —e  ]  +  [e  +2e  ]( 3 6) 

2  4  12  12 


This  expression  is  equivalent  to  equation  (16)  even  though  equation  (35)  is 
not  identical  to  equation  (15).  The  decomposition  characterized  by  (35)  can 
be  found  in  figure  3. 
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3K  4-  2L  +  4(y0  +  y0) 

4(.?  +  s  —  2) 

(4yn  +  fr +  2L)T~ 

4(.?2  +  .v  -  1) 

O  - -> 


I  -K  | 


Figure  3.  Alternate  example  of  transition  to  parallel  representation 


The  application  of  this  result  in  a  six  degree  of  freedom  simulation  is 
illustrated  in  figure  4.  At  any  time  step  t  ,  the  above  analytic  result  can 
propagate  the  output  of  the  transfer  function  to  r  .  The  initial  conditions 
come  from  the  results  of  the  previous  time  step,  the  driving  terms  come  from 
the  aerodynamic  motions  of  the  projectile  (inherently  slower  processes)  which 
are  integrated  numerically  during  the  time  step  interval,  and  possibly  other 
transfer  functions. 


< - A  t  — > 


n  —  1  n  n  +  1  n  +  2 

Figure  4.  Iterative  propagation  of  the  solution 
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When  a  result  is  obtained  for  y  at  t  —  r  ,  it  is  saved  to  become  the  new 
initial  condition  to  be  used  for  once  again  propagating  the  piecewise  analytic 
solution  for  y  to  the  next  time  step  t  Note  that  not  only  must  >'(fn)  be 

saved  to  serve  as  an  initial  condition  for  the  next  time  step,  but  higher  order 
derivatives  as  well.  For  the  example  used  in  this  report,  only 

v</  =  <„)  =  —  (37, 

(It 

is  required.  This  is  easily  obtained  by  differentiating  equation  (16)  or  (36). 


For  this  example,  equations  (4),  (5)  and  (16)  are  used  iteratively  as 
follows: 


Step  A:  L  =  KQ  —  KtQ 

K  is  the  slope  of  the  driving  term  at  t 


o' 


K  is  the  value  of  the  driving  term  at  tQ. 
Initially,  take 


>o  =  0 


>0  =  0 
n  =  0 


Step  B:  Let 


n  =  n  +  1 
t  =  A  t 

Find  y  and  y  at  the  end  of  the  time  step,  t  =  A  t . 


(38) 


(39) 

(40) 

(41) 

(42) 

(43) 


■ 


K^t  K+2L  e  r 
y(Ar)  --  -  -  +  1 4y 0—  4y Q—  AT  +2L 

2  4  12 


+ 


K  +  L+  2y0+y0J  (44) 


y(AO 


K  e 
2  6 


At 


e 

4yQ~4y0~K  +  2L  +  K  +  L  +2y0+y0 


(45) 


Step  C:  Let 


y0  =  >(A  t) 
y 0  =  y(A  t ) 


(46) 
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Evaluate  K  and  L  at  f  (these  depend  ultimately  on  air  frame 
motion  which  is  integrated  numerically  elsewhere  in  the  six 
degree  of  freedom  simulation.) 

Step  D:  Repeat  from  step  B  to  advance  to  the  end  of  the  next  time  step. 


In  summary,  if  the  product  of  several  sequential  transfer  functions  can  be 
recast  into  an  equivalent  network  of  transfer  functions  in  parallel,  difficulties 
arising  from  propagation  of  the  signal  through  a  sequential  network  can  be 
eliminated  even  when  relatively  large  integration  time  steps  are  used.  This 
can  be  done  by  making  a  partial  fraction  decomposition  of  the  Laplace 
operator  representation.  This  technique  produces  a  decomposition  that 
always  leads  to  terms  that  can  be  readily  integrated.  The  use  of  a  partial 
fraction  decomposition  also  facilitates  the  treatment  of  noise.  The  next 
section  contains  a  more  general  and  formal  treatment  of  some  aspects  of  these 
techniques. 
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GENERAL  ANALYSIS 


Partial  Fraction  Expansion 


In  this  section,  a  formal  treatment  is  given  for  carrying  out  the  partial 
fraction  decomposition  using  the  limit  method.  The  function  Y(s)  for  which 
an  inverse  transform  is  desired  typically  appears  as  the  ratio  of  two 
polynomials. 


Vis)  j  =  o 

Yis)  =  -  =  - 

G(s) 

2  V 

k  =  0 


i  i  -  1 

as  +  a.  ,.v  +  ....  +  a.s  +  an 

i  i-l  to 


n  n  -  1 

.r  +  j  -f  —  -h  b  s  +  b 

n  -  1  10 

with  n  >  i  and  b  ~  1 . 


(47) 


Note  that  Y(s)  must  be  a  proper  fraction  (  n  >  i  ).  If  this  is  not  the  case,  an 
improper  fraction  can  always  be  reduced  to  a  proper  fraction  by  long 
division. 

It  is  desired  to  decompose  this  single  complex  transfer  function  into  a 
group  of  simpler  parallel  transfer  functions  that  may  be  more  readily 
handled,  as  explained  in  the  previous  section.  The  technique  to  be  used  to 
achieve  parallel  decomposition  is  called  partial  fraction  expansion.  Suppose 
that  the  denominator  of  equation  (27)  can  be  factored  to  yield  the  form 


V{s)  V(s) 

Y(s)  =  -  =  -  (48) 

G(.v)  is -r1)(.v-r2)(.v-r3)  ....  (s-r^) 

where  the  r,  ,  j  —  1  to  n  are  distinct  roots  of  G(s)  =  0.  (A  discussion  of 
obtaining  ths^  roots  r  will  be  treated  subsequently.)  Then  the  partial  fraction 
expansion  far  distinct  roots  is  given  by 
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Y(s)  = 


+ 


4-  ....  + 


(49) 


C 

n 


s  —  r 


1 


s  —  r 


2 


s  — 


r 

n 


The  expansion  coefficients  C  to  C  are  to  be  determined.  Equating  the  right 
sides  of  equations  (48)  and  (49)  gives 


C 


C, 


+ 


+ 


+ 


c 


s  — 


r 

n 


V(s) 


(.v-r1)(j-r2)(j-r3)  ....  (.v-rj 


V(-0 

G(s) 


(50) 


The  i,h  constant  C.  can  be  determined  by  multiplying  both  sides  of 
equation  (50)  by  the  denominator  of  the  ith  fraction  in  equation  (50),  viz. 
(s  —  r.).  Therefore,  the  i*  factor  will  be  canceled  on  the  left  side  of  equation 
(50).  ‘Finally,  if  r.  is  not  a  repeated  root,  take  the  limit  as  s  goes  to  r.  and 
solve  for  C.  .  This  process  is  expressed  mathematically  by  the  following 
relationship. 

V(s)  is  -  r.) 

C  ]jm  (51) 

Gis) 

for  1  =  1,2,3, . ,n. 


Note  that  G(s)  always  contains  a  term  of  the  form  ( s  —  r .)  which  must  be 
removed  by  division  before  going  to  the  limit.  (For  this  reason,  that  root 
cannot  be  a  repeated  root.  See  below  if  this  is  ppt  the  case.)  Since  the 
inverse  Laplace  transform  of  A  /  (j  —  r.)  is  A  e  ,  it  is  always  a  trivial 
matter  to  invert  Y(s)  in  equation  (27)  to  find  y(t), 

yit)  =  L'1  [t  (.v )  J 


17 


-1 


1  2 
+ 

(•v-z-j)  (s-r2) 


+ 


C 

n 

+ 

(,-rn) 


(52) 


where 


=  2  v .  ( t ) 
i=  l 


v.(r)  =  C  L 

i  i 


r  t 

C  Ae  ' 

i 


(53) 


If  Y(s)  contains  repeated  roots  in  the  denominator,  the  procedure  for 
finding  the  C  for  the  partial  fraction  expansion  must  be  modified.  Consider, 

V(s)  V(s ) 

Y(.v)  =  =  (54) 

G(0  (s-r  )m(s-r  Ms-r  )  ....  (s-r  ) 

i  m  i  m  —  m  -*■  /i 


The  root  r  is  repeated  here  m  times,  and  there  are  n  nonrepeated 
(distinct)  roots  making  a  total  of  (m  +  n)  factors  of  G(s).  Expanding  in 
partial  fraction  form,  equation  (54)  becomes 


Y(s)  = 


Vis) 
O  (s) 
C, 


C 


(•v~r1) 


+ 


is-rj 


m  -  1 


4- 


m—  2 


c 


m  -  1 


+ 


+ 


+ 


+  _ + 


(•V-Tj) 

c 


m  —  j  ■+■  1 


(■v-rp 


+ 


(55) 


(^-r  „) 


The  constants  ...  for  the  repeated  roots  are  determined  from  the 
following  relationships: 
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(56) 


'V(s)  , 

Ci=lim  (s~ri) 

s  -  r ,  Cj  {S  ) 


d  r  V(S)  "  ■ 

=  lim  [j-rJ 

s  -r ,  ds  I  6(.0 


1  d  f  V(s)  f  ‘ 

=  lim  l9_rij  (58) 

.i-r,(k-  1)!  jsk  1  1  (7  (.v) 

Expression  (58)  may  be  derived  as  follows:  Consider  equation  (55), 
which  is  rewritten  here  in  more  compact  form  for  convenience. 

j  m  —  n 

Y(s)  =  ^  Cl(s-r1)‘~(m'l)  +  £  C;(.9-r;)-1  (59) 

/  =  1  /  =  m  -  1 

The  terms  in  braces  and  in  the  second  summation  represent  the  distinct  roots 
of  G(s).  It  is  desired^to  obtain  an  expression  for  C  .  To  obtain  this,  multiply 
both  sides  by  (,r  —  r  )  or  1 


V(s) 


G(s) 


/  =  m  -  1 


If  this  is  differentiated,  the  leading  term  C  vanishes.  The  C2  vanishes  when 
the  second  derivative  is  taken,  ai^d  so  forth.  It  is  clear  that  taking  k-1 
derivatives  reduces  the  typical  (j  )  term  to  a  constant  (all  prior  terms 
vanishing),  with  each  differentiation  bringing  down  an  exponent  of  (s  —  r  ) 
that  leads  to  a  factorial  coefficient,  viz. 


d -1  v(i)  l  m  d~l  ..  m+"  d~l 

{s-r.)  =  SC.  ('-rj  +  S  C,  (s-rfis-r,) 

ds  G  (J )  .  ,  ds  .  .  ds 

1*1  /  ■  m  +  1 


(/-l)!  ,_4  _  m  \ 

(*-l)!Ct+  S  C, - (i-r,)  +  2  C, - 

,  ('—*)!  _  («"*  +  !)! 


(j-Tj)  (i-Tj) 
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similar  terms  in  (s—r,) 


(61) 


This  is  almost  what  is  required,  viz.  an  expression  for  C  .  If  at  this  point  the 
limit  is  taken  as  s  goes  to  r.  ,  all  extra  terms  vanish  leaving  the  expression  in 
equation  ( 68) . 


Finding  the  roots 


One  major  obstacle  that  may  arise  in  applying  this  approach  is  the  task 
factoring  the  denominator  G(s),  i.e.,  finding  the  roots  r,,  r.,  — ,  r  for  the 

in  w  1  2  H 

n  order  polynomial  G(s)  in  equation  (27).  In  most  practical  cases,  the 
transfer  function  is  already  in  the  form  of  factors  of  no  higher  than  second 
order.  These  roots  can  be  found  by  using  well-known  formulae  such  as  the 
quadratic  formula,  viz. 


-  b  ±  ^ I?  2  -  4 ac 

s  = 

2  a 


(62) 


where 


1 

as  +  hs  +  c  =  0  f  63) 

Solving  low  order  polynomial  for  roots  is  not  difficult.  However,  if  higher 
orders  of  n  are  encountered,  techniques  of  higher  sophistication  r/pust  be 
used.  There  are  many  published  techniques  which  solve  n  order 
polynomials.  The  “Laguerre's  Method”  for  polynomial  expansion9  was 
adapted  here. 

The  n  roots  r ,,  /• ,  ....,  r  of  the  polynomial  G(s)  are  sought  where 

12  n 

(}(s)  =  (s  -  r  )(s  -  r  )  .  (s  -  r  )  (64) 

12  n 

n  n  -  1  2 

=  .V  +  c  s  +  .  +  C-S  +  c  v  +  r 

n  -  1  2  1  U 

C  does  not  explicitly  appear  since  all  the  coefficients  of  ,v  in  the  factors 
(s  —  r  )  arc  unity.  Two  expressions  will  be  needed  in  the  following 
development.  The  first  is  obtained  by  differentiating  G(s)  with  respect  to  s. 

"  1 

G  (s  )'  =  G(s)  ^  (65) 

s  —  r  . 

7=  1  J 
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rearranging  gives 


i(.v) 


G(s)' 


G(s) 


+ 


+ 


+ 


(66) 


a 


—  r 

i 


where  a  =  s  —  r  (the  interval  separating  the  variable  s  from  the  n  root). 
The  second  expression  is  obtained  by  differentiating  G(s)'  with  respect  to  s 
again.  From  (65) , 


"1  "1 
G(5)"  -  G(sy  2  -  G(s)  2 - 


s  —  r . 
j=l  J 


y-i(5"V 


(67) 


dividing  through  by  G(s)  and  rearranging  gives 


G(.0" 

-  =  J(s)J(s) 

G  (s ) 


n 


S 

i=  i 


2 


(68) 


or 


^  Forman  S.  Acton,  Numerical  Methods  That  Work,  Harper  and  Row,  New  York,  NY,  1970.  Acton  asserts  that 
Laguerre’s  method  always  converges  to  some  root  from  any  starting  approximation,  even  to  a  complex  root  from  a 
real  first  approximation  if  the  implementation  of  the  algorithm  permits  complex  arithmetic. 
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H(s)  =  J(s)  - 


G(s)‘ 


G(s) 


s  ; 

7-l(5"r/ 


2  2 
(•v-r,)  (.v-r,) 


(•v"rn) 


7=1  s  -  r 

j  . 

While  seeking  the  n  root  r  ,  consider  the  standard  deviation  cr  ,  of  all  the 
other  roots. 

"-1  1  1  f""1  1  ]2 

2  7  -  2 
,=1  n~l  *-/=l 

al-\  = - (70) 

n  —  2 

Using  the  definitions  given  in  (66)  and  (69)  to  eliminate  the  summations  in 
(70)  yields 

2 

2  r  1  i  1  f  1 

in- 2)cr  =  |  H(s)  -  -  -  -  J(s)  -  ~  (71) 

n  1  I  o 

L  a2  n  —  1  L  a  * 

Expressed  as  a  quadratic  equation  in  the  reciprocal  for  a,  this  is 

n  1  2  2 
- -  +  2 J(s)—  -  J(s)  +  (n-l)//(.v)  -(«-!)(« -2)  <rn_1  =  0  (72) 


Using  the  quadratic  formula  to  solve  for  1/a  and  then  inverting  gives 


./  ±  ■  (n  -  \)(nH(s)  -  J2)  -  n(n  -  \)(n-2)a2_1  ■ 
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where  the  sign  of  the  denominator  of  this  expression  for  a  is  chosen  to 
minimize  the  absolute  value  of  a.  To  use  (73)  to  solve  for  the  roots 
iteratively,  the  standard  deviation  cr^  is  neglected.  The  justification  for  this 
step  is  discussed  at  the  end  of  this  section.  Dropping  tr^  gives 


a 


n 


1/2 


./  ± 


r 

| (n  ~  1  ){nH  (.?)  -  J2) 


- 


(74) 


where,  from  (66)  and  (69) 

G(J)' 

J(s)  =  (75) 

G(s) 

and 

G(J)" 

H(s)  =  J(s)  -  -  (76) 

G(s) 


The  iterative  algorithm  for  applying  Laguerre’s  method  is  as  follows: 

1.  Select  an  initial  trial  root,  .vQ,  for  equation  (64). 

2.  Use  the  trial  root  to  evaluate  the  expression  for  a,  using  (75)  and 
(76)  in  (74). 

3.  Generate  an  improved  trial  root  by  replacing  .r0  by  .s-Q  -  a. 

4.  Repeat  step  2  to  obtain  an  improved  value  for  a. 

Continue  the  iteration  of  steps  2  and  3  until  convergence  of  a  has  occurred 
(“a”  sufficiently  small.) 


Laguerre’s  method  can  converge  to  a  pair  of  complex  conjugate  roots, 
since  the  radicand  in  equation  (74)  can  be  negative.  Rapid  convergence 
depends  on  a  good  initial  selection  for  a  trial  root.  A  rule  of  thumb  is  to 
select  a  root  that  is  much  larger  than  the  coefficients  in  the  polynomial  G(s). 

At  this  point,  note  that  neglecting  o-  in  equation  (73)  can  be  justified  as 
follows.  If  the  roots  are  equal,  then  cr  is  identically  zero.  This  can  be 
seen  from  equation  (70)  if  the  roots  become  equal  as  a  limiting  process.  In 
the  limit,  all  the  expressions  s  —  r  tend  toward  the  same  small  value,  which 
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will  be  denoted  by  b.  Then  equation  (70)  becomes 


2 


cr 


n  -  1 


n  —  2 


(77) 


or 


2 


O' 


n  -  1 


1 


(n  ~  1) - 


b 


2 


1 


(n-  1) 


n  —  2 


0 


(78) 


On  the  other  hand,  if  at  least  some  of  the  roots  are  distinct,  neglecting 
is  intuitively  equivalent  to  saying  that  the  root  r  is  far  removed  from  the 
other  roots,  which  are  bunched.  Relatively  speaking,  this  is  always  true  as 
the  solution  begins  to  converge. 


A  FORTRAN  program  implementing  this  technique  is  to  be  found  in 
appendix  A. 
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CONCLUSIONS 


The  use  of  piecewise  closed  integration  techniques  can  result  in 
considerable  relaxation  in  integration  time  step  size  and  significant  savings  in 
computer  run  time.  However,  the  user  must  use  an  equivalent  parallel  rather 
than  a  sequential  representation  to  avoid  unphysical  propagation  delay  effects. 
The  se  of  an  equivalent  parallel  network  is  also  convenient  when  treating 
autopilot  noise  sources.  The  partial  fraction  expansion  techniques  described 
here  are  a  convenient  method  for  transforming  a  product  of  transfer  functions 
into  an  equivalent  parallel  network. 
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SYMBOLS 
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VARIABLE 


DESCRIPTION 


Cl,2,3,..,n 

D 

G(s) 

H(s) 

I 

J(s) 

K 

L 

-1 

L 

r  1,2,3,.. 
s 

T 

t 

V(s) 

Y(s) 

y(t) 


Expansion  Coefficients  of  Partial  Fraction  Decomposition 

Coefficient  of  the  Laplace  operator  s  in  the  second-order  lag/ 
oscillator  with  differentiation 

Denominator  of  transfer  function  Y(s) 

G(s)" 

Transfer  function  J(s)  — 

G(s) 

Coefficient  of  the  square  of  the  Laplace  operator  s  in  the  second- 
order  lag /  oscillator  with  differentiation 

G(s)" 

Transfer  function 

G(s) 

Constant  term  in  the  transfer  function  for  the  second-order  lag / 
oscillator  with  differentiation 

Laplace  transform 


Laplace  transform 
Roots  of  G(s) 

Laplace  transform  operator 

Driving  term  in  the  transfer  function  for  the  second-order  lag / 
oscillator  with  differentiation 

Time 

Numerator  of  transfer  function  Y(s) 

V(s) 

Transfer  function,  written  as  ;  Laplace  transform  of  y(t) 

G(s) 

Inverse  Laplace  Transform  of  Y(s) 
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APPENDIX  A 

FINDING  ROOTS  BY  LAGUERRE’S  METHOD 
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uuu  uuuuuu 


PROGRAM  ROOTS 

IMPLICIT  COMPLEX*16  (A-H.O-Z) 

COMPLEX‘8  J 
REAL*8  ACC,RR ,RI 

DIMENSION  C(21),ROOT(20),CD(21),CDD(21),D(21) 

C 

C  THIS  ROUTINE  FINDS  THE  ROOTS  OF  A  NTH  ORDER  POLYNOMIAL  USING 
C  LAGUERRE’S  METHOD.  POINT  OF  CONTACT:  MICHAEL  J.  AMORUSO  OR 
C  ROMEL  CAMPBELL,  ARDEC, DOVER,  NJ  201-724-4523. 

C 

C  ***  DATA  MUST  APPEAR  IN  A  FILE  CALLED  RDATA  IN  FIELDS  20  COLUMNS 
C  ***  WIDE.  SEE  EXAMPLE  OF  FORMAT  BELOW. 

C 

OPEN  (3,FILE=  ’RDATA’) 

REWIND  3 

INITIALIZING  VARIABLES  AND  CONSTANTS  FOR  POLYNOMIAL  EQUATION. 

ZERO  =  DCMPLX(0.,0.) 

DO  5  I  =  1,21 
C(I)  =  ZERO 
CD(I)  =  ZERO 
CDD(I)  =  ZERO 
5  ROOT  (I)  =  ZERO 

TRIAL  ROOT 

S  =  DCMPLX(100.0D0,0.D0) 

CONVERGENCE-TEST  TOLLERANCE 
ACC  =  1.0D-12 

C  EVALUATING  EQUATIONS  FOR  N'TH  ORDER  POLYNOMIAL 
C  AND  ITS  FIRST  AND  SECOND  DERIVATIVES. 

C 

C  ***  THE  ORDER  OF  THE  POLYNOMIAL  "N"  AND  ITS  N+l  COEFFICIENTS  ARE 
C  ***  ARE  READ  IN  HERE.  COFFICIENTS  ARE  READ  FROM  HIGHEST  POWER  OF 
C  ***  VARIABLE  TO  CONSTANT  TERM. 

C 

C  ***  EXAMPLE:  FIND  THE  ROOTS  OF  4TH  ORDER  POLYNOMIAL: 

C 

C  4  3 

C  S+S-2S=0 

c 

C  ***  THE  INPUT  DATA  WOULD  BE  IN  THE  FOLLOWING  FORM: 

C 
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on  nnn  non  nnnnnn 


c 

C  4 

C  1.0 

C  1.0 

C  0.0 

C  1.0 

C  0.0 

c 

READ(3,101)  N 
WRITE(6,103)  N 
DO  9  I=N+1,1,-1 
READ(3,102)  RR,RI 
C(I)  =  DCMPLX(RR.RI) 

IF  (I.NE.l)  WRITE(6,104)  C(I),I-1 
IF  (I.EQ.l)  WRITE(6,107)  C(I) 

9  CONTINUE 
WRITE(6,105) 

NUMBER  OF  ROOTS  TO  BE  FOUND 
K  =  N 


COEFFICIENTS  OF  1ST  DERIVATIVE  OF  POLYNOMIAL 

1  DO  6  1=1,20 
RR=  DFLOAT(I) 

BUF=DCMPLX(RR,ZERO) 

CD(I)  =  C(I+  1)*BUF 

6  CONTINUE 

COEFFICIENTS  OF  2ND  DERIVATIVE  OF  POLYNOMIAL 

DO  7  1=1,20 
RR=  DFLOAT(I) 

BUF=  DCMPLX(RR,ZERO) 

CDD(I)  =  CD(I+1)*BUF 

7  CONTINUE 

EVALUATING  POLYNOMIAL  AND  ITS  FIRST  TWO  DERIVATIVES 

G  =  C(l) 

GD  =  CD(1) 

GDD  =  CDD(l) 

DO  10  I  =  1,K 
G  =  (C(I+1)*(S**I))  +  G 

GD  =  (CD(I+  1)*(S*  *1))  +  GD 
10  GDD  =  (CDD(I+1)*(S**I))  +  GDD 

POLYNOMIAL  ROOT  SOLVING  ROUTINE  USING  LAGUERRE’S  METHOD 
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nnn 


c 

CK  =  DCMPLX  (  DFLOAT(K)  ,  0.D0) 

CK1  =  CK  -  DCMPLX  (  1.0D0  ,  0.D0) 

J  =  GD  /  G 
H  =  J**2  -  (  GDD/G  ) 

DEN  =  CK1*(CK*H-J*J) 

DEN  =  CDSQRT(DEN) 

IF  (CDABS(J  +  DEN).GE.CDABS(J-DEN))  THEN 
A  =  CK  /  (J+DEN) 

ELSE 

A  =  CK  /  (J-DEN) 

END  IF 
S  =  S  -  A 

IF  (CDABS(A/S).LT.ACC)  GOTO  20 
GO  TO  1 
C 

C  REDUCTION  OF  POLYNOMIAL  BY  SYNTHETIC  DIVISION 
C 

C  THE  FACTOR  (X-S)  WILL  BE  DIVIDED  OUT  OF  THE  POLYNOMIAL 
C  TO  REMOVE  THE  ROOT  JUST  OBTAINED.  THEN  THE  REDUCED  POLY- 
C  NOMIAL  CAN  BE  SOLVED  FOR  THE  NEXT  ROOT. 

C  SEE  FORMAN  S.  ACTON,  "NUMERICAL  METHODS  THAT  WORK", 

C  HARPER  AND  ROW,  NEW  YORK,  1970,  PP  181-183. 

C 

20  CONTINUE 
C  RR  =  DREAL(S) 

RR=  S 

RI=  DIMAG(S) 

ROOT(K)  =  S 
WRITE(6,100)  ROOT(K) 

THE  ORDER  OF  THE  REDUCED  POLYNOMIAL  IS  OBTAINED. 

K  =  K-l 

IF  (K.GE.l)  THEN 
DO  30  I=K  + 1,1,-1 
D(I)  =  (D(I+  1)*S)  +  C(I+1) 

30  CONTINUE 
C 

C  THE  D(I)  ARE  NOW  THE  COEFFICIENTS  OF  THE  OLD  POLYNOMIAL 
C  DIVIDED  BY  THE  FACTOR  (X-S),  WHERE  S  IS  THE  ROOT  JUST  FOUND. 
C  THE  C(I)  WILL  BE  SET  TO  THESE  D(I),  WHICH  WILL  BE  NORMALIZED 
C  SO  THAT  THE  HIGHEST  ORDER  TERM  HAS  A  COEFFICIENT  OF  UNITY. 
C 

C  SET  NORMALIZATION  FACTOR. 

C 

Q  =  D(K+  1) 

DO  40  I  =1,K  +  1 
C(I)  =  D(I)/Q 
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40  D(I)  =  0. 

C 

GO  TO  1 
END  IF 
C 

WRITE(6,106) 

STOP 

100  FORMATC  (’.G20.8, \’,G20.8,’)’) 

101  FORMAT(10X,I10) 

102  FORMAT(2D20.0) 

103  FORMAT(///,’  LAUGERRE  METHOD  FOR  SOLVING  FOR  POLYNOMIAL  ROOTS.’, 
*  //,’  THE  ORDER  OF  THE  POLYNOMIAL  IS  ’,13,//,’  THE  TERMS  ARE:’) 

104  FORMAT(’  (’,G20.9,’,’,G20.9,’)  *  X‘*’,I2) 

105  FORMAT(/,’  THE  ROOTS  ARE:’) 

106  FORMAT(//////) 

107  FORMAT(’  (’,G20.9,’,’,G20.9,’)’) 

END 
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APPENDIX  B 

SECOND  ORDER  LAG/OSCILLATOR  WITH  DIFFERENTIATOR 
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In  a  previous  report,1  transfer  functions  were  treated  using  this  analytic 
approach.  Solutions  were  provided  for  lead/lag,  lag,  lag  with  differentiator, 
and  second  order  lag/oscillator.  Sometimes  these  are  sufficient  and  there  is 
no  need  for  parallel  decomposition.  For  completeness,  the  second  order 
lag/oscillator  with  differentiator  is  developed  here.  Consider  the  differential 
equation 


d  v 

dy 

dT 

I 

II 

& 

+ 

i 

1 

Q 

+ 

(B  -1) 

i 

dt  ~ 

dt 

dt 

where  both  /  and  K  are  nonzero.  The  derivative  of  the  driving  term  appears 
on  the  right  side  as  an  input,  which  is  undesirable.  A  standard  technique 
introduces  an  auxiliary  variable,  say  Z,  and  substitutes  a  pair  of  coupled 
differential  equations  in  the  variables  Z  and  y  that  do  not  contain  any 
derivatives  of  the  driving  term  T .  These  equivalent  differential  equations  are 

dZ 

=  Ky  (B  -2) 

dt 

and 

dy 

I  +  Dy  =  T  -  Z  (B  -3) 

dt 

The  equivalence  of  (B  -2)  and  (B  -3)  to  (B  -1)  can  be  readily  verified  by 
differentiating  (B  -3)  with  respect  to  t  and  substituting  (B  -2)  into  the  result. 


To  implement  the  solution  of  the  coupled  equations  (B  -2)  and  (B  -3),  an 
equation  for  Z  in  which  y  has  been  eliminated  is  sought.  This  is  done  by 
using  (B  -3)  to  eliminate  y  in  the  equation  (B  -2).  This  yields 


dZ 

=  Ky 

dt 

\k  1  ! 

dy 

=  |— 1 
Id  J  1 

T  -  Z  -  I  — 

dt  - 

^  Eugene  M.  Friedman  and  Michael  J.  Amoruso,  "An  Analytical  Modularized  Treatment  of  Autopilots  for  Guided 
Projectile  Simulations,"  US  Army  Armament  Research  and  Development  Center,  Technical  Report  ARLCD-TR- 
85025,  Picatinny  Arsenal,  NJ,  August  1985. 
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(B  -8a) 


(B  -8b) 


The  radicand  D  —  AIK  in  (B  -8)  divides  the  solution  into  classes.  If  the 
radicand  is  nonzero,  the  homogeneous  solution  may  be  written 

Mt  Nt 

Z  (t)  -  A  e  +  B  e  (B  -9a) 

ti 

When  the  radicand  doesn't  vanish,  M  and  N  are  distinct  and  the  two 
exponentials  are  linearly  independent  and  (B  -9a)  is  a  complete  general 
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solution  to  the  homogeneous  form  of  (B  -5).  When  the  radicand  is  positive, 
(B  -9a)  represents  a  pure  damped  solution.  If  the  radicand  is  negative,  then 
M  and  N  are  complex  conjugates  of  one  another.  In  this  case,  the  solution  is 
a  damped  oscillation.  It  may  be  desirable  to  show  this  behavior  explicitly. 
The  solution  (B  -9a)  may  be  rewritten 


ZH(t)  =  A 


sin  [a>/  +  <j>] 


(B  -9b) 


where  \  and  gj  are  the  real  and  imaginary  parts  of  the  complex  conjugate  pair 
M  and  N,  and  C  and  cj)  are  the  constants  of  integration  instead  of  A  and  B  . 

If  the  radicand  vanishes,  (B  -9a)  is  not  general  since  it  does  not  contain  two 
distinct,  linearly-independent  solutions.  Subsjyjution  verifies  that  the  second 
independent  solution  can  be  obtained  from  te  where  M  =  -D/21 .  In  this 

case, 


Ml  Ml 

Z  (t)  =  A  e  +  B  t  e 

H 


(B  -9c) 


The  next  step  is  finding  a  particular  integral  of  (B  -5).  The  damped  positive- 
radicand  solution  will  be  treated  first.  The  method  of  variation  of  parameters 
is  applicable.2  Assume  the  following  form  for  the  solution  to  the 
inhomogeneous  equation,  viz.  (B  -5). 


Ml  Nt 

z  (O  =  U^t)  e  +  u2(t)  e 


(B  -10) 


where  m  and  u  are  to  be  determined.  It  can  be  shown  that  the  u  and  u  of 
( B  - 10)  satisfy 


m  j '  ( t)  =  -  e 


Ml 


u2'(t)  =  +e 


■Nt 


K  T(t) 


I  (N  -  M) 


K  7(0 


I  (N  -  M) 


(B  -11) 


(B  -12) 
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where  the  prime  denotes  the  derivative  with  respect  to  t.2  Integrating  u  '  and 
//  '  yields  the  particular  solution  required  in  order  to  generate  the  complete 
general  solution  to  (B  -5).  However,  this  can  not  be  done  in  closed  form  in 
the  general  case  since  the  form  of  T(t)  is  not  known  analytically  a  priori.  In 
order  to  carry  out  these  integrations  explicitly,  some  generality  can  be  given 
up  that  is  in  the  spirit  of  this  approach.  The  simplest  nontrivial  assumption  is 
that  the  time  step  is  small  enough  to  approximate  T(t)  linearly,  viz. 


T  (r)  =  At  +  b  =  T  '(t  —  t)  +  T 

O  O 


(B  -13) 


where  A,b,T  ,  and  T  '  are  constants.  Integrating  u  '  and  u  2'  with  this 
assumption  yiefds  ' 


K  e 


Ml 


llx(t)  = 


r  -4 1 

I  At  +  b  +  —  | 

MI  (N-M)  l  M  - 1 


-  Mi  , 

K  e  \  T  ' 

- |  T(t)  +  - 

MI  (N-M)  L  M 


(B  -14) 


and 


-K  e 


Nt 


u2(t)  = 


NI  ( N-M ) 


At  +  b  + 

N 


-Nt 


-Ke  \  T  '  ] 

- |  T  +  - | 

NI  (N-M)  L  N  J 


(B  -15) 


Putting  (B  -14)  and  (B  -15)  into  (B  -10)  yields 


2 

C.  R.  Wylie,  Jr.,  Advanced  Engineering  Mathematics,  McGraw  Hill,  New  York,  NY,  1960. 

•5 

With  this  assumption,  (B-l)  can  be  solved  directly  without  (B-2)  or  (B-3)  since  the  right  hand  side  of  (B-l)  is 
replaced  by  the  constant  A  =  T  ’  (i.e.,  /  y  +  Dy  +  Ky  =  A  ). 
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z  (f)  = 

p 


K 

f  r. 

T  \  - 

1  1 

1 

\  +  T  ' 

r  1  1 

J 

I  — 

I  (N-M) 

|  ^  M 

1 

1  2  2 

K  M  N 

x  1 


K  \  ( N-M )  (A/2-  A/2) 

-  j  T - +  T  ' - 

I  (N-M)  l 


2  2 
A/  /V 


a: 


INM 


(N  +  M) 


T  +  T  " 


MN 


(B  -16) 


From  (B  -8),  MN  =  K/  1  and  M  +  N  =  -D/7. 

D 

Z  (t)  =  T  -  —  T  '  (B  -17) 

P  a: 


This  solution  can  be  verified  by  substituting  (B  -17)  into  (B  -5),  keeping  in 
mind  that  T '  is  constant  during  the  time  step.  Since  this  constancy  is  the  only 
assumption  used  in  verifying  (B  -17)  by  substitution  into  (B  -5),  this 
particular  solution  Z  holds  for  all  three  cases  even  though  it  was  formally 
derived  for  only  one  case.  This  result  can  be  combined  with  the  results  for 
the  three  homogeneous  solutions  [viz.,  (B  -9a)  through  (B  -9c)]  and  the 
boundary  conditions  applied  to  determine  the  constants  of  integration. 


POSITIVE  RADICAND 


Combining  (B  -9a)  and  (B  -17)  gives  the  complete  solution  for  the  purely 
damped  case 


Mi  Nl 

Z(t)  =  Ae  +  Be  +  T 


D 

T  ' 
K 


(B  -18) 


43 


dZ{t) 


(B  -19) 


,  Mr  Nt 

=  Z'  =  AMe  +  BNe  +  T 


dt 


The  initial  conditions 


z  = 

Z  at  r 

o 

0 

dZ 

Z'  = 

-  at 

t 

o 

dt 

o 

T  = 

T  at  t 

o 

O 

dT 

T  ’ 

=  T  '  or 

O 

dt 

(B  -20) 


determine  the  constants  of  integration  A  and  B  .  Substituting  into  (B  -18)  and 
(B  -19)  yields 


Mio  Nto  D 

Z  =  A  e  +  B  e  +  T  -  — T  ' 

0  a  o 

K 


(B  -21) 


Mio  Nt 

Z'  =  AMe  +  B  N  e  +  T 


(B  -22) 


solving  simultaneously  yields 


A  - 


Z' 

O 

( 

1 

—  ] 

■ 

z  - 

-  T  +  T  '  \ 

—  + 

o 

N 

o  o  ^ 

N 

K  J 

. 

Mt 


1  - 


M 


N 
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'1  D  JL _ 

T  +  T  '  —  +  — 

o  o  xi 

M  K  1  _  _ 


(B  -23) 


Finally,  y(t)  can  be  obtained  by  substituting  (B  -19)  into  (B  -2),  viz., 


^  Mi  Nt 

y(t)  =  ~  AMe  +  BNe  +  T  ’ 
K 


(B  -24) 


NEGATIVE  RADICAND 


Combining  (B  -9b)  and  (B  -17)  give  the  complete  solution  for  the  damped 
oscillatory  case 


Z{t)  =  Ce  sinfcor  +  4>]  +  T  —  T  ' 

K 


(B  -25) 


dZ(t) 


=  Z'  =  Ce  |  X  sin[wt  +  <{>]  +  u>  cos[o>t  +  <j>]  }■  +  T  '  (B  -26) 


D  V4 IK  -  D 

where  X  =  -  and  oo  =  [See  (B  -8)].  See  (B  -20)  for  the 

21  21 

initial  conditions.  The  constants  of  integration  C  and  <j>  can  be  determined 
now  from 


Z  =  Ce  sin[wr  +  4>]  +  T  -  ~T  ' 

o  o  ~  o  o 


(B  -27) 
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(B  -28) 


Xr.  ( 

Z'  =  Ce  -j  X.  siniojr  +  <bj  +  <*>  cos[co/  +  cj>] 

o  o  o 


[  +  T  ' 


From  ( B  -28) , 


[  Z'  -  T  '  ]  e 


-  X/ 


C  = 


\  sin[wr  +  <j>]  +  co  costw/1  +  c|>] 


Substituting  (B  -29)  in  (B  -27)  gives 


D 


<j>  =  tan 


-  1 


to  (Z  -  T 

o  o 


K 


T  '  ) 

O 


D 

LZ'  -  T  '  -  A(Z  ~  T  +  ~ T  '  ) 

o  o  o  o  o 

K 


—  dit 


(B  -29) 


(B  -30) 


Finally,  y(t)  can  be  obtained  by  substituting  (B  -26)  into  (B  -2),  viz., 


y(t)  = 


K 


( 


j 

Ce  |  \  sin [cu r  +  cj>J  +  to  cos[o>f  +  cj>] 


+  T  ' 


(B  -31) 


ZERO  RADICAND 


Here  M  =  —D/21  [see  equation  (B  -8a)].  Combining  (B  -9c)  with  (B  -17) 
gives  the  complete  solution  for  this  case. 


Mi  Mi  D 

Z(t)  =  Ae  +  Bte  -FT—  T  ' 

K 


(B  -32) 
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1 


dZ(t) 


(B  -33) 


Mt  Mt  Mt 

=  Z'  =  AMe  +  BMte  +  Be  +  T  ’ 
dt 

The  constants  of  integration  can  be  determined  now  from  the  initial 
conditions. 


Mt  Mt  u 

Z  =  Ae  +  Bt  e  +  T  -  ~T  '  (B  -34) 

o  o  o  o 

K 

Mt  Mt  Mt 

Z'  =  AMe  °  +  BMt  e  °  +  Be  °  +  T  '  (B  -35) 

o  o  o 

Solving  (B  -34)  and  (B  -35)  simultaneously  yields 


A 


e 


-  Mt 


O 


D  MD 

Z  -T  +  T  '  +  t  [MZ  ~Z'  -MT  + - T  '  +T  '  ] 

O  O  O  O  O  O  0  o  o 

K  K 


(B  -36) 


and 

MD 

MZ  -  Z'  -  MT  +  - T  '  +  T  ' 

o  o  o  o  o 

K 


(B  -37) 


Finally,  y(t)  can  be  obtained  by  substituting  (B  -33)  into  (B  -2),  viz., 

1  f  Mt  Mt  Mt 

y(t)  —  [  AMe  +  BMte  +  Be  +  T 

K 


(B  -38) 


Initially,  Z  can  be  taken  to  be  zero  on  the  very  first  integration  time  step  and 
Z'  can  be  derived  from  the  initial  condition  y  =  y  at  by  means  of 

o  o  0 

equation  ( B  -2) . 
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